Compare Pseudobulk and Cell Level DE Results
Zhaoheng Li, Wei Sun
2025-08-14
library(Seurat)
library(ggpubr)
library(viridis)
library(tidyr)
library(dplyr)
library(tidyr)
library(readxl)
library(stringr)
library(gridExtra)
library(grid)
library(purrr)
library(tidyverse)
out_dir = "/fh/fast/sun_w/zhaoheng_li/MSK_DEC24/Analysis/outs/differential_expression"
celllevel_out_dir = file.path(out_dir,"cell-level")
pseudobulk_out_dir = out_dir
fig_dir = file.path(out_dir,"compare-DE")for(tp in c("D-1","D3","D7","D11","D18")){
print(paste0("time point: ", tp))
pseudobulk_res = read.delim(file.path(pseudobulk_out_dir, paste0("DE_pseudobulk_", tp, ".tsv")),
header = TRUE, sep = "\t", stringsAsFactors = FALSE)
wilcox_res = read.csv(file.path(celllevel_out_dir,paste0("DE_cell-level_",tp,".csv")),row.names = 1)
merged <- pseudobulk_res %>%
select(celltype, genotype, gene, pvalue) %>%
inner_join(
wilcox_res %>% select(celltype, genotype, gene, p_val),
by = c("celltype", "genotype", "gene")
) %>%
mutate(
pvalue = pmax(pvalue, .Machine$double.xmin),
p_val = pmax(p_val, .Machine$double.xmin),
nl10_pb = -log10(pvalue),
nl10_wc = -log10(p_val)
)
for(ct in unique(merged$celltype)){
print(paste0("cell type: ", ct))
fig = merged %>%
filter(celltype == ct) %>%
ggplot(aes(x = nl10_pb, y = nl10_wc)) +
geom_point(alpha = 0.6, size = 1) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
stat_cor(method = "spearman", label.x = 1, label.y = 19, size = 3) +
coord_cartesian(xlim = c(0, 20), ylim = c(0, 20)) +
labs(
x = "-log10(pseudobulk p-val)",
y = "-log10(wilcoxon p-val)",
title = paste(tp, ct, sep = "_")
) +
theme_classic() +
facet_wrap(~genotype)
print(fig)
ggsave(file.path(fig_dir,paste0(tp,"_",ct,".png",sep='')),fig,width = 12,height = 12)
}
}## [1] "time point: D-1"
## [1] "cell type: ESC"
## [1] "time point: D3"
## [1] "cell type: DE"
## [1] "cell type: ESC_DE"
## [1] "time point: D7"
## [1] "cell type: Liver"
## [1] "cell type: PFG"
## [1] "time point: D11"
## [1] "cell type: PP"
## [1] "time point: D18"
## [1] "cell type: Ductal"
## [1] "cell type: EnP"
## [1] "cell type: SC-alpha"
## [1] "cell type: SC-beta"
## [1] "cell type: SC-delta"
## [1] "cell type: SC-EC"
## [1] "cell type: Stromal"
Session information
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 5461632 291.7 11419752 609.9 11419752 609.9
## Vcells 42015379 320.6 108915709 831.0 108915709 831.0sessionInfo()## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.11.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/Los_Angeles
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] lubridate_1.9.3 forcats_1.0.0 readr_2.1.5 tibble_3.2.1
## [5] tidyverse_2.0.0 purrr_1.0.2 gridExtra_2.3 stringr_1.5.1
## [9] readxl_1.4.3 dplyr_1.1.4 tidyr_1.3.1 viridis_0.6.5
## [13] viridisLite_0.4.2 ggpubr_0.6.0 ggplot2_3.5.1 Seurat_5.1.0
## [17] SeuratObject_5.0.2 sp_2.1-4
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.16.0 jsonlite_1.8.8
## [4] magrittr_2.0.3 spatstat.utils_3.0-4 farver_2.1.2
## [7] rmarkdown_2.26 ragg_1.3.1 vctrs_0.6.5
## [10] ROCR_1.0-11 spatstat.explore_3.2-7 rstatix_0.7.2
## [13] htmltools_0.5.8.1 broom_1.0.5 cellranger_1.1.0
## [16] sass_0.4.9 sctransform_0.4.1 parallelly_1.37.1
## [19] KernSmooth_2.23-22 bslib_0.7.0 htmlwidgets_1.6.4
## [22] ica_1.0-3 plyr_1.8.9 plotly_4.10.4
## [25] zoo_1.8-12 cachem_1.0.8 igraph_2.0.3
## [28] mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
## [31] Matrix_1.7-0 R6_2.5.1 fastmap_1.2.0
## [34] fitdistrplus_1.2-1 future_1.33.2 shiny_1.8.1.1
## [37] digest_0.6.35 colorspace_2.1-0 patchwork_1.2.0
## [40] tensor_1.5 RSpectra_0.16-1 irlba_2.3.5.1
## [43] textshaping_0.3.7 labeling_0.4.3 progressr_0.14.0
## [46] timechange_0.3.0 fansi_1.0.6 spatstat.sparse_3.0-3
## [49] httr_1.4.7 polyclip_1.10-6 abind_1.4-5
## [52] compiler_4.4.0 withr_3.0.0 backports_1.4.1
## [55] carData_3.0-5 fastDummies_1.7.3 highr_0.10
## [58] ggsignif_0.6.4 MASS_7.3-60.2 tools_4.4.0
## [61] lmtest_0.9-40 httpuv_1.6.15 future.apply_1.11.2
## [64] goftest_1.2-3 glue_1.7.0 nlme_3.1-164
## [67] promises_1.3.0 Rtsne_0.17 cluster_2.1.6
## [70] reshape2_1.4.4 generics_0.1.3 gtable_0.3.5
## [73] spatstat.data_3.0-4 tzdb_0.4.0 hms_1.1.3
## [76] data.table_1.15.4 car_3.1-2 utf8_1.2.4
## [79] spatstat.geom_3.2-9 RcppAnnoy_0.0.22 ggrepel_0.9.5
## [82] RANN_2.6.1 pillar_1.9.0 spam_2.10-0
## [85] RcppHNSW_0.6.0 later_1.3.2 splines_4.4.0
## [88] lattice_0.22-6 survival_3.6-4 deldir_2.0-4
## [91] tidyselect_1.2.1 miniUI_0.1.1.1 pbapply_1.7-2
## [94] knitr_1.46 scattermore_1.2 xfun_0.44
## [97] matrixStats_1.3.0 stringi_1.8.4 lazyeval_0.2.2
## [100] yaml_2.3.8 evaluate_0.23 codetools_0.2-20
## [103] cli_3.6.2 uwot_0.2.2 systemfonts_1.1.0
## [106] xtable_1.8-4 reticulate_1.36.1 munsell_0.5.1
## [109] jquerylib_0.1.4 Rcpp_1.0.12 globals_0.16.3
## [112] spatstat.random_3.2-3 png_0.1-8 parallel_4.4.0
## [115] dotCall64_1.1-1 listenv_0.9.1 scales_1.3.0
## [118] ggridges_0.5.6 leiden_0.4.3.1 rlang_1.1.3
## [121] cowplot_1.1.3